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Abstract 

Numerous studies have identified large quantum mechanical efi'ects in the dynamics of liquid 
water. In this paper, we suggest that these effects may have been overestimated due to the use of 
rigid water models and flexible models in which the intramolecular interactions were described using 
simple harmonic functions. To demonstrate this, we introduce a new simple point charge model 
for liquid water, q-TIP4P/F, in which the 0-H stretches are described by Morse-type functions. 
We have parameterized this model to give the correct liquid structure, diffusion coefficient, and 
infra-red absorption frequencies in quantum (path integral-based) simulations. The model also 
reproduces the experimental temperature-variation of the liquid density and affords reasonable 
agreement with the experimental melting temperature of hexagonal ice at atmospheric pressure. 
By comparing classical and quantum simulations of the liquid, we find that quantum mechanical 
fluctuations increase the rates of translational diffusion and oricntational relaxation in our model by 
a factor of around 1.15. This effect is much smaller than that observed in all previous simulations of 
simple empirical water models, which have found a quantum effect of at least 1.4 regardless of the 
quantum simulation method or the water model employed. The small quantum effect in our model 
is a result of two competing phenomena. Intermolecular zero point energy and tunneling effects 
destabilize the hydrogen bonding network, leading to a less viscous liquid with a larger diffusion 
coefficient. However this is offset by intramolecular zero point motion, which changes the average 
water monomer geometry resulting in a larger dipole moment, stronger intermolecular interactions, 
and slower diffusion. We end by suggesting, on the basis of simulations of other potential energy 
models, that the small quantum effect we find in the diffusion coefficient is associated with the 
ability of our model to produce a single broad 0-H stretching band in the infra-red absorption 
spectrum. 



I. INTRODUCTION 



The potential influence of quantum mechanical zero point energy and tunneling effects 
on the properties of liquid water has been appreciated since at least the 1970s, ^ and modern 
computational studies of the problem can be traced back to the pioneering path integral 
molecular dynamics and Monte Carlo simulations of the liquid that were performed in the 
mid-1980s. ^'^ However, the effects of quantum mechanical fluctuations on the properties of 
the liquid have still not entirely been resolved, and the importance of liquid water in a 
wide variety of chemical, biological, and geological contexts is continuing to motivate new 
research on this problem.^ 

The early path integral simulations mentioned above demonstrated that quantum me- 
chanical effects weaken the hydrogen-bonding network in liquid water and make it less struc- 
tured than its classical counterpart.^'^ More recently, the centroid molecular dynamics^'^ 
(CMD) and ring polymer molecular dynamics^'^ (RPMD) approximations have been used 
to study the role of quantum mechanical effects in the dynamics of the liquid. These 
simulations have found that the rates of translation diffusion and orientational relaxation in 
ambient water increase by a factor of approximately 1.5 when quantum mechanical fluctua- 
tions are included. Similar quantum effects have also been obtained using other simulation 
methods, including the Feynman-Hibbs (FH) method^^ and the Feynman-Kleinert linearized 
path integral (FK-LPI) approach, and for a variety of different water models. Overall, the 
general concensus that has emerged from these studies, namely that the disruption of the 
hydrogen-bonding network caused by quantum fluctuations leads to a less viscous liquid 
with faster translational diffusion and orientational relaxation, would appear to be quite 
well estabhshed.^"^^ 

There are however some reasons to believe that this may not be the whole story. First, 
most previous studies that have assessed the role of quantum fluctuations in liquid water 
have used empirical potential models that were parameterized on the basis of classical sim- 
ulations. Quantum simulations of such models lead to a "double counting" of the quantum 
effects in the structure and dynamics of the liquid. This typically results in, for example, 
a computed diffusion coefficient for the room temperature liquid that is much larger than 
the experimental value, indicating that such models are unsuitable for quantum simula- 
tions. Second, many of the empirical water models that have been studied using quantum 
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methods assume that the water monomers can be treated as rigid bodies. Although 
this is convenient from the point of view of computational efficiency, it disregards the im- 
portant role played by intramolecular flexibility in the liquid structure, dynamics, ^^"^^ and 
thermodynamics. Finally, of those flexible water models that have been studied using 
quantum dynamics methods, the majority describe the intramolecular flexibility using sim- 
ple harmonic potentials. ^'^^'^^ The significant anharmonicity of the 0-H vibration of the 
water monomer, which is evident from red-shifting and broadening of the 0-H stretching 
band in the infra-red (IR) absorption spectrum of liquid water, is thereby neglected. 

In the present paper, we shall address these issues by introducing a new empirical water 
model that has been specifically parameterized to reproduce the liquid structure, diffusion 
coefficient, and IR absorption frequencies in quantum, rather than classical, simulations. 
Our four-site q-TIP4P/F water model is both fiexible and anharmonic in the intramolecular 
0-H stretching vibration, and we believe that it provides a more realistic description of the 
effective interactions in the bulk liquid than any model that has been parameterized on the 
basis of classical simulations. We shall also compare and contrast our results with those of 
the recently developed three-site q-SPC/Fw quantum water potential, and demonstrate 
that our model gives a better description of several properties including the melting point 
and the temperature dependence of the liquid density. 

More importantly, we shall show that the quantum mechanical effects in our q-TIP4P/F 
water model are much smaller than have been observed in any previous simulation. For 
example, the diffusion coefficient in our quantum liquid is only around 15% larger than it is in 
the corresponding classical liquid, rather than ~ 50% larger as observed in previous quantum 
dynamics simulations of other water models. This finding can be rationalized in terms of 
the competition between intermolecular and intramolecular quantum effects, and it is a 
direct result of using an anharmonic intramolecular potential. We find that intermolecular 
quantum fluctuations speed up the dynamics of the liquid by disrupting the hydrogen- 
bonding, in agreement with previous quantum simulations. However, quantum fluctuations 
in the intramolecular degrees of freedom result in a larger average water monomer dipole 
moment in the quantum liquid than in the classical liquid. This larger dipole moment 
increases the strength of the intermolecular interactions and retards the molecular motion, 
thereby acting in competition with the intermolecular quantum effect. This competition 
between intermolecular and intramolecular quantum effects has not been described before 
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and it is the main new result of the present work. 

In Section II, we briefly describe the quantum water models we have studied, and in 
Section III we outline our simulation procedures. In Section IV we describe the results 
of extensive quantum simulations of the static and dynamic properties of our q-TIP4P/F 
water model, paying specific attention to the comparison with both experimental results 
and the results of simulations with the previously-developed q-SPC/Fw model.^^ In Section 
V, we compare the results of classical and quantum simulations of the q-TIP4P/F model 
and examine the competition between intra- and intermolecular quantum effects. Section 
VI ends the paper with some concluding remarks. 



II. QUANTUM WATER MODELS 

It is well established that when simulating inhomogeneous environments such as the 
surfaces of water and ice one should use a water model that includes an explicit treatment 
of electronic polarization.^^ However, when simulating the bulk liquid one can often get 
away with using a simple point charge model that has been parameterized to capture the 
polarization in a mean-field sense, ^^'^^ and since this is computationally more convenient it 
is the approach we have adopted in the present paper. 

Of the large number of simple point charge models that have been developed for liquid 
water, the one we have chosen as the starting point for the present work is the four-site, 
TIP4P/2005 classical rigid water model of Abascal and Vega.^^ Unhke most other simple 
point charge models, this has the advantage that it gives a reasonable description of the 
ice/water phase diagram in classical Monte Carlo simulations, correctly predicting the rel- 
ative stabilities of a variety of ice polymorphs. This is an appealing feature of the model 
even for liquid simulations, because one would expect the local structures of the various 
ice polymorphs to be visited in accordance with their relative stabilities at 298 K during a 
simulation of the ambient liquid. 

The TIP4P/2005 model is a pairwise additive intermolecular potential for the interaction 
between water molecules of the form^^ 



i .i>i 
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where Vij is the distance between the oxygen atoms and is a distance between partial 
charge sites in molecules i and j. Two positive charges of magnitude Qm/'^ are placed on 
the hydrogen atoms of each molecule and a negative charge of —qu is placed at a point Tm 
a fraction 7 along the vector connecting the oxygen atom to the center-of-mass of the two 
hydrogens: 

rM = 7ro + (l-7)(rHi+rHj/2. (2) 

In order to add intramolecular flexibility to this model, we have used a quartic expansion 
of a Morse potential to describe the stretching of the 0-H bonds and a simple harmonic 
potential in the bond angle, 



Mntra — 



where 



Von{r) = 



7 



(3) 



(4) 



Here rn and are the two 0-H distances and 9i is the H-O-H bond angle in water molecule 
i. The use of Eq. (4) rather than the full Morse potential VoH(r) = Dr[l- e^^'-^'-'^^i)]^ 
avoids dissociation events, which would not be described correctly using such a simple model 
in any case. An anharmonic description of the 0-H bonds is needed to recover the single 
broad absorption band that is seen in the 0-H stretching region of the liquid water IR 
spectrum, whereas previous simulations have shown that a simple harmonic potential suffices 
to capture the H-O-H bending band.^'^^ 

In contrast to this functional form, the three-site q-SPC/Fw quantum water model of 
Voth and co-workers^^ places the negative charge on the oxygen atom (7 = 1) and uses a 
harmonic description of the monomer stretching and bending modes: 

Von{r)^]^K{T-T,^f. (5) 

As has been shown previously, and we shall demonstrate again below, this purely harmonic 
intramolecular potential does not provide a very good description of the 0-H stretching band 
in the vibrational spectrum of liquid water. However, since the q-SPC/Fw model was also 
parameterized on the basis of quantum simulations, it does provide a convenient benchmark 
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against which to compare our present model. 

In total, the interaction potential in Eqs. (1) to (4) contains nine parameters. The in- 
termolecular parameters are £, cr, qu and 7, all of which we have simply fixed at their 
values in the TIP4P/2005 potential, and the intramolecular parameters are Dr, ar, Tcq, 
ke, and ^eq- These intramolecular parameters were optimized in an extensive series of cal- 
culations to give good agreement with the experimental structure, self-diffusion coefficient, 
and vibrational absorption frequencies of the liquid in quantum mechanical (path integral) 
simulations. The final parameters of the resulting q-TIP4P/F model are fisted along side 
those of the q-SPC/Fw modeP=^ in Table I. 

III. QUANTUM SIMULATION METHODS 
A. Path Integral Molecular Dynamics 

We have used the path integral molecular dynamics (PIMD) method^^ to calculate a 
variety of structural and thermodynamic properties of the q-SPC/Fw and q-TIP4P/F mod- 
els, including radial distribution functions, dielectric constants, liquid densities, and melting 
points. Most of these calculations were straightforward and require little additional ex- 
planation. For example, the 0-0, 0-H and H-H radial distribution functions of the two 
liquids were obtained from 250 ps NVT path integral simulations in the presence of an 
Andersen thermostat. Dielectric constants were obtained from much longer (10 ns) NVT 
simulations in order to allow for complete dielectric relaxation, and liquid densities at 1 atm 
pressure were obtained from 10 ns NPT simulations in order to fully converge the average 
over density fiuctuations. These latter simulations were performed in the presence of both 
an Andersen thermostat and an isotropic Berendsen barostat.^^ 

The melting point calculations were a little more complicated and do need some more 
explanation. These were done by performing direct coexistence simulations of the water-ice 
interface under atmospheric pressure. Initial hexagonal ice configurations were generated 
by placing the oxygen atoms at their crystallographic sites. The hydrogen atom positions 
were determined using the Monte Carlo procedure of Buch et al.^'^ in such a way that the 
Bernal-Fowler rules^^'^^ were satisfied and the total dipole moment of the simulation cell 
was exactly zero. The initial ice configuration was then equilibrated in the presence of 



6 



an Andersen thermostat and an anisotropic Berendsen barostat for 50 ps,^^ before placing 
the secondary prismatic (1210) face of the ice cell in contact with an equilibrated water 
simulation. 

In total, the coexistence simulations consisted of 696 water molecules, with 360 initially 
in the ice phase and 336 in the liquid. The combined ice/water system was simulated for 10 
ns in the presence of an Andersen thermostat and an anisotropic Berendsen barostat. Both 
the number density profile of oxygen atoms along the axis perpendicular to the ice/water 
interface and the total potential energy of the system were used as order parameters to mon- 
itor the extent of melting or freezing during the simulation. The simulation was halted if 
complete melting or freezing was observed to have occurred. The melting temperature was 
determined to within 1 K for the q-TlP4P/F model and to within 5 K for the q-SPC/Fw 
model by repeating the whole procedure at different temperatures and using a bisection 
procedure to home in on the melting point. 

B. Ring Polymer Molecular Dynamics 

We have also used the ring polymer molecular dynamics (RPMD) method^'^ to calculate 
several dynamical properties of the room temperature q-TlP4P/F liquid, including the self- 
diffusion coefficient and various orientational relaxation times. The diffusion coefficient was 
obtained from the time integral of the RPMD velocity autocorrelation function, 

1 1"^ 

^ = 3 Cv.v(i) dt, (6) 

and the Z-th order relaxation times for various axes 77 of the water molecule were obtained 
from the time integrals of the corresponding orientational correlation functions, 

/•oo 

t]'^ / cl{t)dt. (7) 

Both of these calculations have been described in detail in recent papers and we have used 
the same procedures in the present study.^^'^^'^^ The velocity autocorrelation function Cv.v(^) 
in Eq. (6) was calculated for 2 ps by time averaging over 100 consecutive NVE RPMD 
trajectories of length 4 ps, with a re-sampling of the momenta from the Maxwell distribution 
between each trajectory. This procedure was repeated twenty times to produce an average 
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value for and a standard error in the diffusion coefficient. Tfie orientational correlation 
functions Cj{t) were calculated up to 10 ps by averaging over 500 NVE trajectories of 
length 20 ps, again with a re-sampling of the momenta between each one. Where necessary, 
an exponential tail was fit to the correlation function beyond 10 ps before evaluating the 
integral in Eq. (7). 

Voth and co-workers have recently calculated both D and the same r/'s as us for their 
q-SPC/Fw water model using the centroid molecular dynamics (CMD) method. This is 
closely related to RPMD and it gives similar results for systems like liquid water that are 
nearly classical. We ourselves prefer to use RPMD to calculate D and r/' because it gives a 
more accurate approximation to both c^.^{t) and cj'(t) in the short time limit, which is the 
only limit in which anything has been proven rigorously about the relative accuracies of the 
two approximations for a general potential.^ 



C. Partially-Adiabatic Centroid Molecular Dynamics 

In a recent paper, we have shown that IR spectra calculated using RPMD exhibit 
spurious peaks arising from the internal modes of the harmonic ring polymer. These peaks 
appear at high frequencies (above 1300 cm^^ in liquid water at 298 K), and so do not interfere 
with the calculation of diffusion coefficients and orientational relaxation times, both of which 
are zero-frequency spectral components [see Eqs. (6) and (7)]. However, the spurious peaks 
do contaminate the 0-H stretching band of the vibrational spectrum of liquid water, which 
is centered at around 3400 cm~^. 

A simple fix to this problem is to adjust the elements of the Parrinello-Rahman^^ mass 
matrix so as to shift the spurious oscillations beyond the spectral range of interest, as is 
done in the partially-adiabatic centroid molecular dynamics (PA-CMD) method of Hone, 
Rossky and Voth.^^ We have therefore used this method in the present work to calculate the 
vibrational spectra of the q-TIP4P/F and q-SPC/Fw water models. As we have recently 
argued in more detail,^^ it is convenient to choose the elements of the Parrinello- Rahman 
mass matrix so that the internal modes of the ring polymer are shifted to a frequency of 

„n/(n-l) 
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where j3 — l/k^T and n is the number of beads in the ring polymer {n — 32 in the 
calculations reported below). This prescription gives Q/27rc = 7414 cm~^ at 298 K, which is 
sufficiently far beyond the spectral range of interest in the liquid water vibrational spectrum 
that increasing Vt any further does not have any significant impact on the results. "^^ It also 
ensures that the PA-CMD simulations can be performed with the same integration time step 
as is used in RPMD.^^ 

The IR absorption spectra of the q-TIP4P/F and q-SPC/Fw models were calculated 
using the formula^^'^*^'^^ 



Here n{u}) is the frequency-dependent refractive index of the liquid, a{uj) is a Beer-Lambert 
absorption coefficient, V is the volume of the simulation cell, c is the speed of light in a 
vacuum, and Cij,.^{t) is a canonical (Kubo-transformed^^'^^) dipole autocorrelation function. 
This was calculated up to 10 ps by time-averaging along 300 NVE PA-CMD trajectories 
of length 20 ps, with a re-sampling of the momenta between each trajectory. The Fourier 
transform in Eq. (10) was then performed numerically. Further details of this approach to 
calculating the vibrational spectrum of liquid water are given in ref. 38. 

D. Additional Computational Details 

Unless stated otherwise, all of our liquid simulations were performed at a temperature of 
298 K and density of 0.997 g cm~^ with 216 water molecules in a cubic simulation box. Pe- 
riodic boundary conditions were applied using the minimum image convention. Short-range 
interactions were truncated at 9 A and Ewald summation was employed to calculate the 
long-range electrostatic interactions.^^ 32 ring polymer beads were used in all simulations. 
A ring polymer contraction scheme with a cut-off value of cr = 5 A was employed to reduce 
the electrostatic potential energy and force evaluations to single Ewald sum, thereby signif- 
icantly speeding up the calculations.^^ Several classical simulations were also performed for 
comparison by collapsing the ring polymer to a single bead (n = 1). The ring polymer time 
evolution was performed analytically in the normal mode representation using a multiple 




(9) 



where 




(10) 
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time-step algorithm in which the intermolecular forces were updated every 0.5 fs and the 
intramolecular forces every 0.1 fs. In all simulations, the system was equilibrated for 100 ps 
in the presence of an Andersen thermostat^^ before the accumulation of any averages. 



IV. VALIDATION OF THE Q-TIP4P/F MODEL 

We shall now validate the q-TIP4P/F water model by comparing its predictions for 
various properties of liquid water with experimental measurements when the properties 
are calculated using the path integral methods described above. We shall also compare the 
results of this model with those of the recently-developed q-SPC/Fw quantum water model.^^ 
The competition between intra- and intermolecular quantum effects in the dynamics of the 
q-TIP4P /F model is a separate issue that we shall investigate by making a comparison with 
classical molecular dynamics simulations in Section V. 



A. Static equilibrium properties 

Let us begin by considering the ensemble-averaged water monomer properties of the q- 
TIP4P/F and q-SPC/Fw water models given in Table II. It is clear from this table that 
there are significant differences between the average quantum geometries of the q-TIP4P/F 
and q-SPC/Fw water monomers in the room temperature liquid. The bond length is greater 
in q-SPC/Fw by about 0.04 A, and the bond angle is greater by about 1.5°. Despite the 
smaller partial charges on the q-SPC/Fw hydrogen atoms, these geometric changes result 
in an average monomer dipole moment in q-SPC/Fw water which is greater than that in 
q-TIP4P/F water by about 5%. 

Table II also reports average values of the tetrahedral quadrupole moment 

Qt = 9i^^9yy^ (11) 

where Qxx and Qyy are the quadrupole tensor components along the axis joining the hydrogen 
atoms and normal to the molecular plane. Qt gives an indication of the size of quadrupolar 
interactions in the liquid and has recently been shown to correlate strongly with the classical 
melting points of rigid-body water models.^^ The value of (Qt) for q-SPC/Fw is about 20% 
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smaller than that for q-TIP4P /F and this has significant implications for the melting points 
of the two models as we shall see below. 

Figure 1 compares the 0-0 radial distribution functions (RDFs) of the two models with 
the results of the most recent neutron scattering experiment of Soper.^^ Both calculated 0-0 
RDFs possess too much structure in the first peak, indicating that short-range interactions 
between nearest-neighbour oxygen atoms are slightly too strong; however, both models are 
consistent with the experiment to within its reported error bars.^^ The slightly increased 
structure observed in the first peak for q-TIP4P/F compared to q-SPC/Fw is as a result of 
its larger tetrahedral quadrupole moment (see Table II) . To improve the level of agreement 
in this feature, a more sophisticated treatment of the short-range interactions would be 
necessary; recent work using Gaussian charge sites and Buckingham potentials has indeed 
been shown to improve matters, albeit with an increase in the number of parameters 
defining the interaction potential. Beyond the first maximum in goo{f): the q-TIP4P/F 
model does somewhat better than the q-SPC/Fw model in reproducing the oscillations in 
the experimental RDF. 

Our computed 0-H and H-H RDFs for the two water models are compared with the 
neutron scattering results in Figs. 2 and 3. One sees from these figures that the experimental 
intramolecular 0-H and H-H peak positions are correctly reproduced by the q-TIP4P/F 
model. This is clearly a result of our having included these RDFs in the parameterization of 
the model. In contrast, the q-SPC/Fw model predicts a somewhat larger 0-H bond length 
and a larger bond angle than is seen in the experimental data, leading to intramolecular 0-H 
and H-H peaks at radii that are shghtly too large. This is consistent with the ensemble- 
averaged monomer geometries of the two water models given in Table II. 

Constant-pressure PIMD simulations of both quantum water models give calculated den- 
sities at 298 K and atmospheric pressure that are in good agreement with experiment. Our 
calculated density for the q-SPC/Fw model is also in good agreement with that of an ear- 
lier path integral simulation. However, as shown in Fig. 4, the density of the q-SPC/Fw 
water model is only really satisfactory at 298 K, the temperature at which the model was 
parameterized.^'^ As the temperature is increased or decreased, the agreement with exper- 
iment deteriorates. In particular, our calculations indicate that the density maximum for 
the q-SPC/Fw model occurs at 235 K, around 40 K lower than the experimental value, and 
at a density which is 2% greater than the experimental density maximum. 
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In contrast, we find that the q-TIP4P/F model gives a very good description of the 
temperature-dependence of the hquid density. Throughout the entire range between the 
melting and vapourization points of liquid water at 1 atm, the difference between the q- 
TIP4P/F and experimental densities is less than 0.3%. This ability of the q-TlP4P/F 
model to reproduce the experimental density curve is clearly inherited from the properties 
of the underlying rigid-body TIP4P/2005 model.^^ One reason why it is desirable to get the 
density curve right before simulating the ambient liquid is that it has been suggested to have 
some bearing on the ability of a water model to provide a good description of hydrophobic 
solvation. ^'^ 

The melting temperature of the q-TlP4P/F model obtained from our coexistence simu- 
lations of the ice-water interface was found to be 251 ± 1 K, around 22 K lower than the 
experimental value. In contrast, our coexistence simulations of the q-SPC/Fw model were 
found to give a melting point of 195 ± 5 K, about 78 K lower than experiment. The larger 
error bars associated with this value are a result of the slower motion at such low tempera- 
tures, which increases the time scale required to observe melting or ice lattice formation at 
the interface. 

To put these numbers in context, classical simulations of common rigid water models 
have been found to give melting points which range from about 146 K for TIP3P to 272 K 
for TIP4P/Ice (which was specifically parameterized to reproduce this property). Classical 
simulations have also found a strong correlation between the melting point and the tetra- 
hedral quadrupole moment Qt, with a larger quadrupole moment increasing the relative 
stability of the ice phase leading to a higher melting point. '^^ The higher quantum melting 
point of the q-TlP4P/F model is therefore associated with its larger average tetrahedral 
quadrupole moment (Qt)- It is also clearly inherited from the properties of the underlying 
rigid-body TIP4P/2005 model, which melts in classical simulations at 252 K.^^ Three-site 
models from the SPC and TIP3P famihes generally have smaller tetrahedral quadrupole 
moments and lower melting points,^^'^^ and our computed quantum melting point for the 
q-SPC/Fw model is entirely consistent with this. 

Previous classical simulations have also suggested that it is not possible to reproduce the 
experimental difference of 3.98 K between the temperature of maximum density and the 
melting point of hexagonal ice using simple point charge models. The present quantum 
mechanical simulations of the q-TIP4P/F model give a 28 K difference between these two 
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temperatures, which is in the middle of the 21-37 K range of differences found in classical 
simulations of similar models. So the inclusion of quantum mechanical effects does not 
seem to help matters. It appears that an explicit treatment of electronic polarization will be 
needed to reproduce the small temperature difference between the temperature of maximum 
density and the melting point of water. This is consistent with the finding that the average 
dipole moment of the molecules in ice is significantly larger than that of the molecules in 
liquid water, indicating that significant charge redistribution occurs upon freezing. ^^'^^ 

Another property that is likely to benefit from an explicit treatment of electronic polar- 
ization is the dielectric constant. The final row of Table II gives the dielectric constants 
obtained in the present path integral simulations of the q-TIP4P/F and q-SPC/Fw models 
at 298 K and 0.997 g cm~^. These were calculated as 

^r-'^ + ^{{l^-l^)-{l^)-{l^)), (12) 

where /u is the total dipole moment of the system. Our result for the q-SPC/Fw model 
is consistent with the value obtained by Paesani et al}^ in an independent path integral 
simulation (e^ = 86 ± 4). Both models give only qualitatively reasonable values for the 
dielectric constant: q-TIP4P/F underestimates the experimental value by 22% while q- 
SPC/Fw overestimates it by 15%. These errors are not surprising given the lack of explicit 
electronic polarization in either model. Note also that the larger dielectric constant of the 
q-SPC/Fw model is consistent with its larger average monomer dipole moment. 



B. Dynamical properties 

Turning now to dynamical properties, Table III gives the results of our RPMD simulations 
of the diffusion coefficient and various orientational relaxation times for the q-TIP4P/F 
model at 298 K and 0.997 g cm~^. The experimental diffusion coefficient of liquid H2O was 
employed as a target during the parameterization of the model, and our calculated value of 
0.221 ps~^ is seen to be in good agreement with the experimental value^^ of 0.230 
ps~^. Perhaps more importantly, since it provides a direct experimental indication of the 
role of quantum mechanical effects, we obtain equally good agreement with the experimental 
H2O/D2O diffusion coefficient ratio. ^® 
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One should however bear in mind that diffusion coefficients are quite sensitive to finite 
size effects/^'^^'^^'^^ and the simulations reported in Table III were performed for a system 
of only 216 molecules. We have therefore performed some additional RPMD simulations 

of larger systems (containing 343 and 512 H2O molecules), and extrapolated the diffusion 
coefficient to the limit of infinite system size using the formula^^'^^ 

where 77 is the shear viscosity, L is the length of the simulation box, and ^ is a numerical 
constant that depends on the geometry of the simulation = 2.837 for a cubic box). This 
was found to give D{oo) = 0.245 ± 0.007 ps-^ for the q-TIP4P/F model, within 7% of 
the experimental result. For comparison with this, the diffusion coefficient of the q-SPC/Fw 
model calculated using CMD has been reported to be 0.24 ps~^ for a system of 216 
water molecules. We would expect that this value would also increase somewhat when 
extrapolated to the limit of infinite system size. 

The remaining rows of Table III give the orientational relaxation times for three different 
axes of the water molecule (the H-H axis, one of the two equivalent 0-H axes, and the 
dipole axis) obtained from RPMD simulations of the q-TIP4P/F model, along with the 
CMD results for the q-SPC/Fw model from Paesani et al}'^ A previous study has shown 
that these orientational relaxation times are less sensitive to system size effects than the 
diffusion coefficient and so in this case a simulation of a system of 216 water molecules is 
perfectly satisfactory.^^ Comparing the calculated relaxation times with the values available 
from NMR and IR relaxation experiments, we find that both the q-TIP4P/F and q- 
SPC/Fw models give a reasonable description of orientational relaxation in bulk liquid water 
in quantum simulations of the room temperature liquid. 

Figure 5 compares our calculated PA-CMD dipole absorption spectra of the q-TIP4P/F 
and q-SPC/Fw water models with the experimental spectrum of Bertie and Lan.^^ Our spec- 
trum for the q-SPC/Fw model is in good agreement with the fuUy-adiabatic CMD spectrum 
reported previously by Paesani et al}^ The two calculated spectra clearly reproduce the gen- 
eral features of the experimental spectrum, with 0-H stretching absorptions above ~3000 
cm~^, a water bending band at around ~1600 cm^^, and intermolecular librational features 
below ~1000 cm~^. However, the peak at ~200 cm~^ is absent from both of the simulated 
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spectra. This peak arises from the low-frequency modulation of dipole-induced dipole inter- 
actions which are clearly not present in simple point charge models. The relative intensities 
of the intramolecular bending and stretching bands also disagree with experiment in both 
simulations, and this again arises from the neglect of electronic polarization. The intensity 
of an absorption band is proportional to the square of the derivative of the dipole mo- 
ment along the corresponding vibrational coordinate, and neglecting the (quite significant) 
induced dipole contribution to this gives the wrong intensity. 

Putting these obvious deficiencies of both models aside, the q-TIP4P/F model clearly 
does a better job of reproducing the experimental vibrational spectrum of liquid water than 
the q-SPC/Fw model. The position and width of the broad 0-H stretching baud in the 
experimental spectrum are both well reproduced by the q-TIP4P/F simulation, whereas 
the q-SPC/Fw model predicts distinct antisymmetric and symmetric stretching peaks that 
are blue-shifted relative to the experimental absorption maximum by around 80 and 190 
cm~^. The appearance of these separate peaks is a direct consequence of the use of a 
harmonic 0-H stretching potential in the q-SPC/Fw model [see Eq. (5)]. The position 
and width of the experimental intramolecular bending band are also fairly well reproduced 
by the q-TIP4P/F simulation, although in this case there is a slight red shift relative to 
the experimental absorption maximum of some 40 cm~^. The red shift in the case of the 
q-SPC/Fw simulation is around 190 cm~^. 

All of these observations are clearly consistent with the fact that the experimental absorp- 
tion frequencies of liquid water were included in the data set that was used to parameterize 
the q-TIP4P/F model. As a further check on the reliability of the model for describing 
the vibrations in the liquid, we have therefore performed some additional PA-CMD simula- 
tions of heavy water, for which the model was not explicitly parameterized. The resulting 
q-TIP4P/F and q-SPC/Fw spectra are compared with the experimental IR spectrum^^ of 
D2O in Fig. 6. One sees from this figure that the q-TIP4P/F model again does a better job 
of reproducing the positions and widths of the bands in the experimental spectrum than the 
q-SPC/Fw model, and indeed that the agreement between the q-TIP4P/F and experimental 
spectra is no worse in this case than it was for H2O. 
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C. Summary 



Combining all of the above results, it is clear that the q-TIP4P/F model does a good job 
of reproducing a wide variety of static and dynamic properties of liquid water in quantum 
mechanical (path integral) simulations. For all of the properties we have considered, the 
results of the q-TIP4P/F model are at least as good as those of the recently developed 
q-SPC/Fw quantum water model, and for some properties (such as the melting point, the 
temperature dependence of the liquid density, and the vibrational frequencies of the liquid) 
they are significantly better. This is not especially surprising, as these properties were either 
included in the parameterization of the q-TIP4P/F model or inherited from the behavior of 
the underlying TIP4P/2005 potential. What is surprising, however, is just how small the 
quantum mechanical effects in the diffusion and orientational dynamics of the liquid turn 
out to be when simulated with this new quantum water model, as we shall now discuss. 

V. COMPETING QUANTUM EFFECTS 
A. Computational results 

The only obvious way to estimate the effect of quantum mechanical fluctuations on the 
dynamics of liquid water is to compare classical (1 bead) with path integral {n bead) simula- 
tions. The results of such a comparison for the diffusion coefficient of the q-TIP4P/F hquid 
are shown in Table IV along with those of similar comparisons that have been made in the 
past for a variety of other water models. ^^^^'^"^'^^ For the q-TIP4P/F model, we find that the 
magnitude of the quantum effect, defined as Dqm/-Dci, is 1.11. This is significantly smaller 
than the effect of quantum fluctuations on the diffusion coefficient that has been found pre- 
viously for any other water model, and indeed the average value of Dqm/-Dci obtained from 
the earher simulations listed in Table IV is 1.46. 

A similar result is found for orientational relaxation. Table V compares the orienta- 
tional relaxation times for three different axes of the water molecule obtained from quantum 
(RPMD) and classical simulations of the room temperature q-TIP4P/F liquid. The quantum 
effect, deflned as Td/r^rn, is seen to be around 1.18 for all three flrst order relaxation times 
and around 1.25 for the second order relaxation times. These effects are again much smaller 
than those seen in all previous studies of orientational relaxation in liquid water, 
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which have typically found the quantum effect Tci/Tqm to be around 1.5. Clearly, then, the 
effect of quantum fluctuations on both the translational diffusion and the orientational re- 
laxation of our q-TIP4P/F model is significantly smaller than has been seen previously for 

any other water model. But why? 

B. Analysis and discussion 

A clear clue as to what is going on is provided by the ensemble-averaged water monomer 
properties of the q-TIP4P/F model obtained in classical and quantum (PIMD) simulations 
of the room temperature liquid (see Table VI) . The average 0-H bond length of the water 
molecules in the quantum liquid is greater than in the classical liquid, and the average bond 
angle is slightly smaller. As a result of these geometric changes, the average water molecule 
dipole moment in the quantum liquid is larger than that in the classical liquid by about 
1.6%. 

A similar effect has been seen previously in an ab initio PIMD study of liquid water, 
which found the average water dipole moment to be larger than that obtained in a classical 
simulation by as much as 4%. The increase in the 0-H bond length in the quantum simu- 
lation is also consistent with a recent experimental investigation which has shown that the 
average 0-H bond length in liquid H2O is larger than the average 0-D bond length in liquid 
D20.^^ The results for the q-TIP4P/F model in Table VI are not therefore unreasonable, 
and indeed they arise from an intuitively rather obvious effect: zero point fiuctuations in 
the anharmonic 0-H stretching coordinate increase the average 0-H bond length and give 
the water molecule a larger average dipole moment than it would otherwise possess. 

The implications of all this for the magnitude of the quantum effect in the diffusion and 
orientational relaxation of the liquid are as follows. A larger average dipole moment leads to 
stronger intermolecular interactions which slow down the translational and orientational mo- 
tion of the molecules in the liquid. However, there is also a competing effect: intermolecular 
quantum fiuctuations disrupt the hydrogen bonding network leading to more rapid diffusion 
and orientational relaxation. For a water model in which this competing effect is the only 
one in operation, such as a rigid-body model or a flexible simple point charge model with a 
harmonic 0-H stretching potential, Dqm/-Dci and Td/rqin arc typically found to be around 
1.5. But for the q-TIP4P/F model, in which both effects contribute and oppose one another. 
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Dqm/Dci and t^Jt^^ are both reduced to around 1.15. 

A simple way to verify this explanation is to turn off the intramolecular flexibility in the 
q-TIP4P/F model and see what happens to the ratio of the quantum and classical diffu- 
sion coefficients. This can be done by performing simulations of the liquid using the same 
intermolecular interactions as in the q-TIP4P/F potential, but with a fixed intramolecu- 
lar geometry. We have therefore calculated both quantum (RPMD) and classical diffusion 
coefficients in this way, using the SHAKE^^ and RATTLE^^ algorithms to constrain the 
intramolecular geometry of the q-TIP4P/F model to its equilibrium value. This was found 
to give = 0.440(2) ps"^ and = 0.308(3) ps"^ in simulations of 216 wa- 
ter molecules at room temperature, and thus a quantum effect in the diffusion coefficient 
of Dfpa/Dc\ = 1.43. When the competition from intramolecular quantum fluctuations is 
removed, we therefore recover the well-established result^^ that intermolecular quantum 
fluctuations speed up the dynamics of the liquid by a factor of at least 1.4. 

The only remaining puzzle in Table IV concerns the results for the TTM2.1-F model,™ 
which has been found to give a quantum effect in the diffusion coefficient of D^^/ — 1.50 
in CMD simulations of a system of 216 water molecules. This is puzzling because TTM2.1- 
F is a ffcxiblc and polarizable Thole-type^^~^^ model with an anharmonic intramolecular 
potential.^^ Why, then, is there not a competition between intramolecular and intermolecular 
quantum effects in this model that brings Dqm/ D^i closer to one? 

We beheve that this is because the TTM2.1-F monomers dissociate to the correct products 
in the gas phase (uncharged radicals). Consequently, the molecular dipole moment in this 
model only shows a slight increase as the 0-H bond length is increased from its equihbrium 
value, and indeed less so than in a typical simple point charge model [see, in particular. 
Fig. 2(c) of ref. 75]. As a result, the energetic cost of stretching an 0-H bond in the liquid 
is not compensated by such a favorable increase in intermolecular interactions as it is in our 
q-TIP4P/F model, and the 0-H bonds in the TTM2.1-F hquid remain fairly close to their 
equilibrium values.™ The monomers in this liquid therefore only access regions close to the 
minimum of the intramolecular potential energy surface where the harmonic character of 
the potential dominates. This is confirmed by a recent CMD simulation of the IR spectrum 
of the TTM2.1-F model, which shows two distinct absorption bands in the 0-H stretching 
region (a feature that is typical of harmonic models). Thus the TTM2.1-F model behaves 
more like the harmonic q-SPC/Fw model than our anharmonic q-TIP4P/F model in terms 
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of both its diffusion coefficient and its IR spectrum. 

A modification to tfie TTM2.1-F model fias recently been developed with the specific 
aim of giving a better description of the vibrational spectra of water clusters and the bulk 
liquid/^ This new TTM3-F model has a molecular dipole moment that increases more rapidly 
as the 0-H bond is stretched, and it therefore behaves more like our q-TlP4P/F model. In 
particular, it gives a single broad band in the 0-H stretching region of the IR spectrum 
in both classical and PA-CMD simulations,^^ and more recent RPMD simulations^^ have 
shown that it gives a quantum effect in the diffusion coefficient of D^^/Dd ~ 0.9. The fact 
that this ratio is actually ^655 than one implies that there must be a very strong competition 
indeed between intramolecular and intermolecular quantum effects in the TTM3-F model, 
with the intermolecular quantum effect dictating the final result. 

More generally, it seems from all of the simulations we have seen that water models which 
give rise to two distinct peaks in the 0-H stretching region of the IR absorption spectrum 
also give a fairly large quantum effect in the diffusion coefficient (-Dqm/-Dci > 1-4), whereas 
models which give rise to a single broad absorption band in the 0-H stretching region give 
a much smaller quantum effect (Dqin/-Dci < 1-2). Since the experimental IR spectrum of 
liquid water shows a single broad band in the 0-H stretching region, this suggests that 
the net quantum effect in the diffusion (and also the orientational relaxation) of the room 
temperature liquid may not be that large after all. 

VI. CONCLUDING REMARKS 

In this paper, we have developed a new potential energy model for liquid water, q- 
TIP4P/F, which gives good agreement with a wide variety of static and dynamic properties 
of the liquid in quantum mechanical (path integral) simulations. We have also analysed the 
effect of quantum mechanical fluctuations on the dynamics of this model in some detail and 
identified two competing contributions to the quantum diffusion coefficient. Intramolecular 
zero point fluctuations increase the average 0-H bond length and the average molecular 
dipole moment, leading to stronger intermolecular interactions and slower diffusion, while 
intermolecular quantum fluctuations disrupt the hydrogen bonding network leading to more 
rapid diffusion. In our q-TIP4P/F model, these two effects nearly cancel one another, leading 
to a comparatively small net quantum effect on the diffusion coefficient. We have argued 
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that the same is hkely to be true for other flexible water models that give a single broad 
absorption band in the 0-H stretching region of the liquid water IR spectrum, as is seen 
experimentally. 

Finally, although we have focussed on the effect of quantum fluctuations on the dynamics 
of the liquid, we believe that the influence of these fluctuations on various static equilibrium 
properties is also likely to be less pronounced for the q-TIP4P/F model than for most other 
water models. To provide just one example of this, we have calculated the classical melting 
points of the q-TIP4P/F and q-SPC/Fw models by repeating the ice- water coexistence 
simulations described in Sec. III. A with the ring polymers replaced by classical particles. 
For the q-SPC/Fw model, ^'^ this was found to give a classical melting point of 222 ± 2 K, 
some 27 K higher than the quantum melting point. For the q-TIP4P/F model, the classical 
melting point was found to be 259 ± 1 K, just 8 K above the quantum result. This is clearly 
a smaller quantum effect, and it is also more consistent with the experimentally-observed 
difference of 4 K between the melting points of H2O and D2O. 
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q-TIP4P/F q-SPC/Fw 



e I kcal mol-i 0.1852 0.1554 

G I k 3.1589 3.1655 

gM / |e| 1.1128 0.84 

7 0.73612 1.00 

Dr I kcal mol"^ 116.09 

ar I 2.287 

kr I kcal - 1059.162 

req / A 0.9419 1.0000 

ke I kcal mor^ rad-^ 87.85 75.90 

^eq / degrees 107.4 112.0 

TABLE I: Parameters in the q-TIP4P/F and q-SPC/Fw (ref. 13) quantum water models. 



25 



q-TIP4P/F q-SPC/Fw Experiment 



iron) I A 
(6'hoh) / degrees 
(/i) /D 
(Qt) / D A 



0.978(1) 
104.7(1) 
2.348(1) 
2.403(1) 



1.019(1 
106.2(1 
2.465(1 
2.009(1 



105.1'^ 



0.97' 



■a 



p29s / g cm 
Ptmd / g cm 

^TMD / K 
^melt / K 



-3 



-3 



0.998(2) 
1.001(2) 
279(2) 
251(1) 
60(3) 



1.002(2 
1.020(2 
235(2) 
195(5) 
90(3) 



0.997^ 
1.000^ 
277.13^ 
273.15'' 



78.4'^ 



TABLE II: Static equilibrium properties of the q-TIP4P/F and q-SPC/Fw quantum water models 
obtained from PIMD simulations. The standard errors in the final digits are given in parentheses. 
[ji) and {Qt) are the ensemble-averaged values of the molecular dipole and tetrahedral quadrupole 
moments, p298 is the liquid density at 298 K and 1 atm pressure, Ptmd is the density at the 
temperature of maximum density, Tmeit is the melting point, and is the dielectric constant. The 
experimental values are from (a) ref. 48, (b) ref. 77, (c) ref. 78. 
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q-TIP4P/F 


q-SPC/Fw 


Experiment 


L>H20 / ps^i 


0.221(1) 


0.24(1) 


0.230'^ 


Dd^o / A2 ps-i 


0.172(1) 


- 


0.177'^ 


-DH20/-DD2O 


1.28 


_ 


1.30« 


Ti"^ / ps 


5.40(4) 


3.8(1) 


_ 


/ ps 


5.10(4) 


4.1(1) 




/ ps 


4.64(3) 


4.7(1) 




^2™ / ps 


2.22(2) 


1.85(5) 


1.6-2.5^ 


^2°" / ps 


1.90(2) 


1.70(5) 


1.95'' 


T2'^ / ps 


1.52(1) 


1.60(5) 


1.90^ 



TABLE III: Dynamical properties of liquid water at 298 K and 0.997 g cm~^ (and heavy water 
at 1.107 g cm~^) obtained from RPMD simulations of the q-TIP4P/F model. D is the diffusion 

coefficient and is the Z-th order orientational relaxation time for molecular axis r/. The CMD 
results for the q-SPC/Fw model from ref. 13 are provided for comparison. The experimental values 
are from (a) refs. 55 and 56, (b) refs. 59-62. 
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Potential Method Dg / ps~^ / ps~^ Dq_m/Dd 



q-TIP4P/F 216 RPMD 

SPC/E 216 RPMD" 

SPC/F 216 RPMD^ 

SPC/F 125 CMD"^ 

TTM2.1-F 216 CMD^ 



0.199(2) 
0.242(3) 
0.279(2) 
0.30(2) 
0.150(5) 
0.358 



0.221(1) 
0.343(2) 
0.400(3) 

0.42 
0.225(5) 

0.548 



1.43 



1.42 



1.11 



1.40 



1.50 



TIP4P 



256 



CMD^ 



1.53 



TABLE IV: Classical versus quantum diffusion coefficients of several liquid water models, including 
the q-TIP4P/F model developed in this work. The second column gives the number of water 
molecules used in each simulation, and the third indicates the approximate quantum dynamical 
method employed; further details can be found in (a) ref. 11, (b) ref. 37, (c) ref. 9, (d) ref. 14 and 
(e) ref. 10. Where given, the number in parentheses indicates the standard error in the final digit. 
The last column gives the magnitude of the quantum effect, defined as Dqm/^ci- The average 
quantum effect from the earlier studies in the table is 1.46. 
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Tel / pS 


^qm / PS 


'^cl/''"qm 


T-HH 


6.3(1) 


5.40(4) 


1.17 


T-OH 


6.0(1) 


5.10(4) 


1.18 




5.48(7) 


4.64(3) 


1.18 


'2 


2.72(4) 


2.22(2) 


1.23 


T-OH 


2.38(4) 


1.90(2) 


1.25 


^2 


1.92(3) 


1.52(1) 


1.26 



TABLE V: Classical versus quantum (RPMD) orientational relaxation times for the q-TIP4P/F 

water model at 298 K and 0.997 g cm~'^. The numbers in parentheses are the standard errors in 
the final digits. The final column gives the ratio of the classical and quantum relaxation times as 
a measure of the magnitude of the quantum effect in orientational relaxation. 
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Classical 


PIMD 


(roH) / A 


0.963(1) 


0.978(1) 


(6'hoh) / degrees 


104.8(1) 


104.7(1) 




2.311(1) 


2.348(1) 



TABLE VI: Water monomer properties for the q-TIP4P/F model in classical and quantum (PIMD) 
simulations. Calculations were performed at a temperature of 298 K and density of 0.997 g cm~^. 
The standard errors in the final digits are given in parentheses. 



30 



o 

D5 







' 1 

- 


' 1 


1 1 1 1 1 

- 

— Experiment 


l.\ 




-- q-TIP4P/F 


nv 

rv 

f V 




n-SPC/Fw 


H V 

- h Vl 

H V 

P V 












H 

A 1 


1 


I.I. 



7 



r/A 



FIG. 1: (Color online.) Oxygen-oxygen radial distribution functions of the q-TIP4P/F and q- 
SPC/Fw quantum water models obtained from PIMD simulations at 298 K and 0.997 g cm~^. 
The experimental radial distribution function from ref. 48 is shown for comparison (along with its 
associated error bars). 
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FIG. 2: (Color online.) Oxygen- hydrogen radial distribution functions of the q-TIP4P/F and q- 
SPC/Fw quantum water models obtained from PIMD simulations at 298 K and 0.997 g cm~^. 
The experimental radial distribution function from ref. 48 is shown for comparison (along with its 
associated error bars). The inset shows the intramolecular 0-H peak at distances close to 1 A. 
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FIG. 3: (Color online.) Hydrogen- hydrogen radial distribution functions of the q-TIP4P/F and 
q-SPC/Fw quantum water models obtained from PIMD simulations at 298 K and 0.997 g cm~^. 
The experimental radial distribution function from ref. 48 is shown for comparison (along with its 
associated error bars). 



33 



1.04 



1.02 - 



E 
o 



0.98 - 



0.96 - 



0.9^ 





1 1 1 

— Experiment 




.. • o q-TIP4P/F 




• q-SPC/Fw 






- 






•_ 

1.1,1 



100 



250 



300 



350 



T/K 



FIG. 4: (Color online.) Densities at 1 atm pressure of the q-TIP4P/F and q-SPC/Fw water 
models obtained from PIMD simulations. The experimental density curve from ref. 77 is shown 
for comparison. The red and blue curves through the calculated points are simply guides for the 
eye. 
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FIG. 5: IR absorption spectra of the q-TIP4P/F and q-SPC/Fw water models obtained from PA- 
CMD simulations at 298 K and 0.997 g cm~'^. The upper panel shows the experimental spectrum 
from ref. 63. 



35 



Experiment 



q-TIP4P/F 



20000 
15000 
10000 
5000 


o 15000 
1 0000 
'S 5000 


150001- q-SPC/Fw 

1 0000 - 
5000 - 











1000 



2000 



(0 / cm 



1 



3000 



FIG. 6: As in Fig. 5, but for D2O at a temperature of 298 K and a density of 1.107 g cm ^. In 
this case the experimental spectrum is only available above 700 cm~^ (ref. 65). 
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